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Abstract 

A fully relativistic quark model is constructed and applied to the 
study of wave-functions as well as the spectrum of heavy-light mesons. 
The free parameters of the model are a constituent quark mass and 
(on the lattice) an adjustable r-parameter in the fermionic kinetic en- 
ergy, while the confinement is introduced via potentials measured by 
MonteCarlo. The results are compared to Monte Carlo energies and 
Coulomb-gauge wave functions. They are in very good agreement with 
the data. A comparison with previous models suggests that we are see- 
ing in the Monte Carlo data the quantum-relativistic derealization of 
the quark due to Zitterbewegung. 
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1 Introduction 



Recent studies in quenched lattice QCD have led to a considerable ad- 
vance in our understanding of meson wave functions - in particular, of the 
relation between the Bethe-Salpeter wavefunction of a heavy-light meson in 
Coulomb gauge QCD and the wavefunctions obtained from a spinless rela- 
tivistic quark model (SRQM) defined by a Hamiltonian of the form |3| 

H x = ^fp 2 +m 2 + V{r) (1) 

where m is a constituent quark mass, and V(r) the confining potential (de- 
termined by Monte Carlo measurements of Wilson line correlations of static 
color sources). 

Wave functions obtained from (|l|) have proved to be enormously useful in 
constructing appropriately smeared lattice operators for heavy-light mesons 
[1], leading to accurate lattice calculations of B-meson properties. They have 
also been recently applied to the extraction of the Isgur-Wise function [|J]. 
Relativistic potential models have also been used to estimate pseudoscalar 
meson decay constants [|| 

Despite the fact that SRQM wavefunctions give a vastly better fit than 
nonrelativistic ones to the meson wavefunctions measured in Monte Carlo 
calculations, some persistent discrepancies in simultaneously describing the 
asymptotic (large distance) behavior as well as the wavefunction at the origin 
suggest that the model defined by Eq([l|) is not capturing all of the essential 
physics, even at the level of a valence quark description. Recall that the 
SRQM of ([l]) has only a single free parameter, the constituent quark mass m, 
as the potential V(r) is determined by Monte Carlo measurements for each 
lattice studied. These discrepancies are not very important in constructing 
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smeared operators for the ground state meson in each angular momentum 
channel, but become very troublesome when one tries to extract excited 
state properties using the multistate formalism of Ref [|J , where admixtures 
of the ground state should be kept to a minimum. 

Our objective in this paper is not only to construct an improved version 
of the SRQM which does a better job in fitting the global behavior of meson 
wavefunctions for different angular momenta and for small as well as large 
distance, but also to provide a clear explanation of the approximations being 
done and the relation of the resulting model with a hypothetical full QCD 
solution of the problem. 

The resulting model extracts, we believe, the full content of the phys- 
ical picture provided by the valence quark description and consistent with 
QCD. The accurate predictions for the wave functions as compared to Monte 
Carlo simulations (see Section 3.1) indicates that Heavy-light mesons can 
be represented reasonably well in terms of this picture. 

The two main effects which emerge from the more complete treatment 
given in Sections 2 and 3 below of the lattice QCD Coulomb gauge Hamil- 
tonian, and which are found to improve considerably the agreement of the 
model with the measured Monte Carlo wavefunctions are 

1. A renormalization of the Wilson r-parameter away from the bare value 
(r=l) used in the Monte Carlo simulations. The sign of this lattice 
effect can be understood already from the one-loop seagull correction 
(see Section 2.1), although the magnitude (as in the case of the quark 
mass correction renormalizing K c ) seems to involve a large nonpertur- 
bative piece. This is reasonable, since a renormalization of r is an ef- 
fect involving all momenta, in particular low momenta where we know 
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perturbation theory fails. Also , one must keep in mind that a one- 
loop calculation in the 4-dimensional Euclidean theory (with at ^ 0), 
will not necessarily give the correct quantitative shift of the spatial 
r-parameter in the Hamiltonian formulation (where a continuum limit 
at — > has implicitly been taken). 

This effect, which should became irrelevant in the continuum limit, 
plays however an important quantitative role improving the agreement 
between model and data for the lattice sizes tested so far (see section 
2.2). 

2. Some of the observed discrepancies between model H\ and the Monte 
Carlo simulations persist, even after the corrections implied in point 
1. These remaining discrepancies are considerably reduced when the 
correct relativistic treatment of the heavy-light system is performed. 
A detailed analysis of the differences between this correct treatment 
and the previous models give rise to a beautiful explanation of this 
new corrections. They turn out to be due to the derealization of the 
light quark known as Zitterbewegung, that, as is well known, arise 
from the inability to localize a relativistic particle in a local unitary 
theory. To my knowledge, these effects are seen for the first time in 
Monte Carlo measured wave functions. 

In section 2, we construct a model that correctly takes into account the 
Wilson lattice fermionic kinetic energy || used in the Monte Carlo simu- 
lations. This model however does not represent an improvement over H\. 
The reason for that is analyzed and as a result a new model arises, incorpo- 
rating the renormalization of the Wilson r-parameter, that does represent 
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an improvement over ([[]). In section 2.2 we compare this new model and 
Hi with the Monte Carlo data. In section 3 we carry out a fully relativistic 
treatment of the problem. In section 3.1 this model is compared with the 
Monte Carlo data. In section 4 we compare the physical content of the three 
models and interpret the differences. In section 5 we present the conclusions 
and discuss upcoming studies. 

2 Improved Treatment of Kinetic Terms 

As was shown in Refs. [[j], ^],the Hamiltonian given by equation (|l]) describes 
very well the results of Monte Carlo calculations of the Coulomb gauge wave 
functions of a heavy-light meson in quenched approximation. In addition to 
practical implications for lattice studies, this model provides a surprisingly 
simple physical picture for the heavy-light mesons, namely, the heavy quark 
acting as a source of the confining Coulomb potential and the light quark 
moving relativistically in this confining field (the relativistic nature of the 
kinetic energy was essential Q in reproducing the large distance behavior 
of the wave function). The real gluons are completely decoupled from the 
quarks except for their role renormalizing the mass. 

In this paper, we will carry the physical picture implied by a valence 
quark model to its limits. The resulting model highly improves the one 
given by Eq- (Jlj) both conceptually and in its predictive power while keeping 
the underlying simplicity. 
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2.1 Using the Wilson Action 

A first, perhaps obvious modification to Hi amounts to replacing the ki- 
netic energy by the lattice Wilson dispersion relation || taking correctly 
into account the specific lattice formulation employed in the simulations. It 
is important in assessing the quantitative validity of the relativistic quark 
model that systematic effects due to lattice discretization be dealt with con- 
sistently both in the model and in the Monte Carlo simulations so that 
deviations between the two may be properly attributed to important phys- 
ical effects rather than lattice artifacts which will eventually disappear in 
the continuum limit. The Monte Carlo calculations Q that constitute the 
'experimental' data were done with a Wilson r parameter equal to one. So 
our new Hamiltonian becomes: 



H' 

where 
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M 2 (q)+£Q? + V(r) (2) 
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M(q) = m+]T(l-cos % ) (3) 
fc=i 

Qk(q) = sinq k (4) 

Although this model is closer to lattice QCD since it contains the cor- 
rect dispersion relation, the corresponding wave functions do not represent 
an improvement with respect to model (|T|). Actually, they magnify the dis- 
crepancies between model Hi and Monte Carlo data. This is at first sight 
very surprising because, as already said, Eq(j2|) is closer to lattice QCD in 
its treatment of the fermionic kinematics than H%. 

The solution to this puzzle comes from a detailed analysis of the renor- 
malization of the parameters of the theory on the lattice. 



More specifically, consider the one loop contribution to the quark self 
energy. On the lattice we have two graphs rather than one (as a consequence 
of the compact representation of the gauge field): 
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(b) 



Figure 1: One loop graphs contributing to the quark self-energy 
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where p = 2 sin ^ and p^ = sinp^J] 

Graph (b) also appears in the continuum while graph (a) is present only 
on the lattice in a compact formulation of the gauge theory. It is precisely 
graph (a) that will provide in the cleanest way the solution to our puzzle, 
as its contribution to the self energy in Coulomb gauge is: 
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where f2 = L 4 , greek indices run from 1 to 4 and roman indices from 1 to 
3 (this convention applies in all equations in this paper) . Eq(^)contains 
1 We are using here the notation of RefQ 



the contribution from the Coulombic instantaneous interaction while Eq(^) 
contains the contributions from the real gluons. 

Writing p 2 as Ylu=i 2(1 — cosp M ), the inverse free propagator becomes 

4 

Ap 1 = m + Ar — r cosp^ + ■ p (8) 

and we immediately realize that the part of Xp 1 ^ proportional to the identity 
matrix (in the Dirac indices) explicitly renormalizes the Wilson r parameter. 
Specifically: 
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For our lattice size, the perturbative r renormalization due to graph (a) are, 
in Coulomb gauge: 
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We shall be comparing RQM models with MonteCarlo data generated 
on a 12 3 x24 lattice at =5.7, corresponding to a naive bare lattice coupling 
$q ~ 1.05. The hopping parameter was k = 0.168. Nonperturbative effects 



may partially be included by using instead the tadpole- improved [ 10 1 defi- 
nition of the coupling, which gives for the (3 value considered a value closer 
to 2.9 for g 2 @. This is the value used in Eqs(0H). 

In our Hamiltonian models we consider of course only r space . This value, 
as we will see in the next subsection correctly predict the sign of the change in 
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r although the magnitude seems to have big nonperturbative contributions. 
Graph (b) also contributes effectively to the r renormalization, but not in 
an explicit way as in the case of the first one. However, in this case the 
numerical contribution is much smaller (as in the case of the mass shift). 

Of course the mass is also renormalized as is well known, and also by 
an amount which is quite a bit larger than the perturbative one-loop shift 
(even with tadpole improved couplings). 

The important point of this calculation is to realize that not only the 
mass but also the Wilson r parameter should be considered as free parame- 
ters, since both of them are dynamically modified, in a nonperturbative way. 

Including this effect, the model acquires the same form as in Eq(0) 



H 2 

but with 



\ 



Af 2 (q) + EQ? + V(r) (13) 

i = l 



M(q) =m + r ^(1 -cos q k ) (14) 
k=l 

We have now therefore two adjustable parameters, m and r. This new 
model, with correctly chosen values for the parameters, represents a sub- 
stantial quantitative improvement over model (Q) as will be shown in the 
next section. We also understand now why Eq.(^) actually works worse 
than Eq(|l|), as H\ is effectively close (in the sense that the fermionic kinetic 
dispersion relation is close to the bosonic one over most lattice momenta) 
to one particular case of the model H2. In fact, it corresponds, for fixed m, 
to r ~ 0.85 as can be seen simply by plotting the corresponding dispersion 
relations. This value, although not optimal, is closer to the optimal choice 
for model ( |l~3|) (see Section 2.2) than the naive unrenormalized choice r = 1. 
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The improvement obtained with Eq (|I3|), although very significant from 
a quantitative point of view for the lattice sizes tested so far , should never- 
theless become irrelevant in the continuum limit, although it is certainly 
relevant in providing accurately smeared meson operators for multistate 
MonteCarlo studies ffl. 

In any case, we have now not only a better model but one that has 
a closer connection to QCD since it contains the actual dynamical QCD 
fermionic kinetic energy . We shall see in the next section that the modi- 
fication in the dispersion formula greatly improves the fit to the measured 
wavefunctions at shorter distance (and in particular at the origin) once the 
m and r parameters are chosen to optimize the fit at medium and large 
distances. 

A fuller description, starting with the Bethe-Salpeter equation (which 
for a light quark propagating in the color field of a static source reduces to 
a Dirac equation) will lead in Section (3) to a model giving similar wave- 
functions, agreeing even more closely with the measured ones. Such a model 
represents a valence quark description of the heavy-light meson that is as 
close to QCD as possible without leaving the physical picture outlined in 
the introduction. 

2.2 Quantitative Consequences of the Improved Potential 
Model 

In order to actually solve for the wave functions of the model, we used the 
same method as in Refs Q We briefly explain it here for completeness. 

The procedure used in a multistate smearing calculation of heavy-light 
meson properties || for generating lattice smearing functions from the RQM 
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is as follows. One obtains orthonormal lattice wavefunctions, which are 
eigenstates of a lattice RQM Hamiltonian defined on a L 3 lattice (with f, f* 
lattice sites): 

Hfftr = Kffti + V{f)8pp (15) 

The eigenstates in a channel of given orbital quantum numbers (S,P,D 
etc) are obtained by applying the resolvent operator E \ H to a source wave- 
functions of the same orbital symmetry. The model at this stage is spinless 
(the measured wavefunctions represent spin-averages of the top two Dirac 
components of the light quark field) so issues of spin-orbit coupling do not 
yet arise (they will be dealt with properly in the full Dirac formalism of 
Section 3). 

In the resolvent approach, S-states are generated by applying the resol- 
vent kernel to a monopole localized at the origin, P-states with a source 
dipole, and so on. At each trial value of the energy E, the norm of the 
resulting state E ^_ H ^^ is evaluated. Obviously 

r' 

when E — > eigenvalue of H. Typically, wavefunctions accurate to 4-5 sig- 
nificant figures are obtained by stopping once this norm exceeds 3000. At 
this point a smearing eigenstate \l4mear(r) is extracted by renormalizing the 
vector E ^_ H ^^ to unit norm. The inversion of E — H is performed by the 
conjugate gradient algorithm, with the multiplication of the kinetic term 
done in momentum space using a fast Fourier transform. 

In the following figures we present the results of our new model as com- 
pared with the old one. As was already mentioned in the previous section, 
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the MonteCarlo data presented in the foolowing figures was generated on a 
12 3 x24 lattice at (3 =5.7, and the hopping parameter was k = 0.168. 

In FigQ we compare the ground state wave functions for the V = 12 3 
case. The values chosen for the constituent mass and r parameters are 
chosen to maximize the agreement (in a mean square sense) between data 
and the respective models in the ground state. As it turned out, the optimum 
constituent masses are very similar to one another and the wave functions 
very insensitive to small changes around the optimum value. We present 
here the results for the same values of the constituent mass. This choice, 
while essentially identical to the optimum cases, helps to appreciate the 
effect of the r renormalization. The case of H2 with r = 1 is also included to 
emphasize the effect of the r renormalization. As we can see, the agreement 
with the Monte Carlo data was already very good for Hi and is further 
improved, specially at the origin by H2. 

But the most important reason for which model H2 was introduced, 
was to capture the lattice artifacts unavoidably present in the Monte Carlo 
data. Only after these artifacts are well under our control can we hope to 
find some physics in the data beyond the one provided by flj. In this sense 
the improvement at the origin is due to the r renormalization as can be seen 
by comparing with the unrenormalized case denoted H2, m = 0.23, r = 1.0 
Also we present in Fig|| a detail of FigQ corresponding to the region of 
distances between R = 1.4 and R = 2.4. Specifically, as can be seen in 
FigQ, at points corresponding to distances R\ = \/3 ( this corresponds to 
the lattice points x\ = ±li ± lj ± lk ) and R2 = 2 (corresponding to the 
point X2 = ±2i ± Oj ± OA; and the points generated by cyclic permutation of 
the coordinates), there is a pronounced 'discontinuity' in the Monte Carlo 
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Figure 2: IS state, L = 12. We see at the origin the improvement of H2 
over Hi when r is renormalized. With r = 1 however, model #2 does a 
poor job showing the necessity of r renormalization. The Monte Carlo wave 
functions were extracted at different time slices. Although all time slices 
gave very similar results, the wave function extracted at the fourth one, 
that we present here, was the one with the best signal to noise ratio. That 
is the meaning of the t4 in the Monte Carlo data point label. 
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data that should of course disappear in the continuum limit. On a finite 
lattice and still not very close to the continuum this discontinuity is easy to 
understand qualitatively: it is due to the fact that under these conditions 
the system responds more naturally in terms of a metric notion of distance 
between two points on a lattice given by some function of the number of 
links between these points (notice that x\ is at 3 links away from the origin 
while #2 is only at 2, in contradistinction with their euclidean distance). In 
figure H we see how model H\ completely ignores this lattice artifact, H2 
with r = 1 is slightly closer, while H2 with r = 0.54 follows almost perfectly 
the discontinuity. 

In Fig[Q] we can better appreciate the large distance region. 

In Fig||], we show the same information as in Fig[Q] but in logarithmic 
scale to appreciate the asymptotic region. As we see, both, H\ and H2 with 
renormalized r do a very good job in this region. 

So, as we have seen, as far as the ground state is concerned, H2 not only 
shows an improvement over H2 specially visible at the origin, but it also 
proved capable of capturing very pronounced lattice artifacts. Both effects 
clearly show the relevance of taking into account the r renormalization. 

Once the values of m and r are specified to reproduce as accurately as 
possible the ground state, we compare now the results for the IP state. In 
this case, we divide the respective wave functions by cos to show only the 
radial dependence. As we can see in Fig||, the model H2 does again a better 
job than H\, although there is still room to improve. The case of H2 with 
r = 1 is not shown since it was in the previous figures only to see the effect 
of renormalizing r. In any case, it again performs worse than H\. 

So we conclude that, although the modifications leading to H2 are only 
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Figure 3: Detail of figure 2. We see here the 'discontinuity' between points 
at distances R\ = 1.73 and R2 = 2. While model H\ completely ignores it, 
and model H2 with unrenormalized r can do just slightly better, model H2 
with the renormalized r almost perfectly follows the discontinuity. 
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Figure 4: Large R region of figure 2. The case H2 with r = 1.0 is not 
displayed to clarify the relevant information. We see that both H\ and H2 
with renormalized r fall very close to the data in this region (notice the 
scale) . 
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Figure 5: IS state L = 12 , logarithmic scale. Both, H\ and with 
renormalized r do a very good job at large distance. 
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Figure 6: IP state , L = 12. The values of the parameters were fixed 
to reproduce as accurately as possible the ground state. We can see the 
improvement of H2 over Hi, but still we have plenty of room to improve. 
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due to lattice artifacts, the quantitative improvement is significant, so the 
value of H2 resides in the fact that it captures a very important lattice 
discretization effect. Nevertheless the improved model is still conceptually 
and quantitatively inadequate. The conceptual inadequacy stems from the 
fact that the relation between the eigenstates of H2 and the spin-averaged 
Bethe-Salpeter wavefunctions in Coulomb gauge is unclear (e.g. the poten- 
tial model ignores antiquarks whereas there are coupled upper and lower 
components in a Dirac formalism). Quantitatively, we shall see that use of 
a full Dirac formalism which is closely related to the Bethe-Salpeter wave- 
function also further improves the agreement with the Monte Carlo results. 
In this full formalism, it will still be important however to include the r- 
renormalization discussed above. 

3 Full Bethe-Salpeter treatment of Heavy-Light 
Wavef unct ions 

As we have seen, the agreement between the wavefunctions derived from 
the Hamiltonian H2 and the Monte Carlo data is quite remarkable; how- 
ever, not only is there still room for further quantitative improvement but 
from a conceptual point of view the connection between these simple mod- 
els and a full hypothetical QCD solution of the meson Coulomb gauge wave 
functions is not completely clear. In another words, it would be nice to have 
a model that works as well as the previous one and in which the nature of 
the approximations being done is completely transparent. In this subsec- 
tion we will construct this model and as a bonus the resulting one will show 
an additional quantitative improvement over H2 with a very nice physical 



19 



interpretation. 

We shall assume that: 

(a) Transverse gluon interactions with the quarks act primarily to renor- 
malize the mass and r parameters in the quark kinetic term. Fock states 
involving real gluons in addition to the valence quarks are neglected. 

(b) The net effect of Coulomb gluon exchange between the light and static 
quarks can be expressed by the potential acting between two infinitely heavy 
color sources. 

More qualitatively, the picture in the back of our mind, supported by the 
comparison with data as will be seen in Section (3.1), consists of the light 
quark moving fast enough for relativistic effects to be important, but on the 
other hand not so fast that the interaction with the static quark cannot be 
accurately described by the energy which would obtain if the light quark were 
held fixed. Alternatively, one might assume that the time scales over which 
the string connecting the quark to the static source responds to changes in 
the light quark position are small compared with the time scales relevant 
for the light quark motion. 

Before we proceed with the derivation of our new model, it will be useful 
to present a brief description of what was actually measured in the Monte 
Carlo simulations of Ref|lJ that constitutes our data. Even though this work 
used a sophisticated multistate smearing method, for our purposes it suffices 
to know that the basic information was extracted from the measurement, in 
quenched lattice QCD, of the Green function: 

Ftf,x,t) = (0\QH(d,th 5 q H (x',t)q H (x,0) l5 Q H (d,0)\0} (17) 

in the limit were the b-quark is taken to be infinitely massive. In this 
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limit, the heavy quark propagator is simply proportional to 1 ~ 2 7 , therefore 
F becomes proportional to the average of the upper two components of 
the light quark propagator in the presence of a color source. From the 
calculation of this object, using the above mentioned multistate smearing 
method (i.e. smearing the source point x of the light quark with Ansatz 
meson wavefunctions derived from Hi) the upper two components of the 
meson wave function were extracted and spin averaged. The result of this 
operation constitutes the data against which we compare our models. 

Taking this into account we will now construct a model that represents 
as closely as possible the quantities measured in the Monte Carlo simulations 
realizing at the same time the physical ideas presented above. 

In a full QCD treatement of the problem at hand, the relevant Bethe- 
Salpeter wavefunction would be 

X (x,t) = (0\q H (x,t)Q H (6,t)\P) (18) 

where |0) is the vacuum, \P) is the meson state (in the center of mass 
frame with energy H\P) = Ebs\P}), and qH,Qn are the light and heavy 
Heisenberg fields. 

In the infinitely massive heavy quark limit, but otherwise still in full 
QCD, Eq.(p^) is best written as, 

X (x,t) = (0\q H (x,t)\P*) (19) 

where |P*) = Qjf(0, t)\P ). This notation emphasizes the fact that in the 
above limit, the heavy quark field is not dynamical. 



As we see, if we were able to calculate exactly Eq.(19) in the context 



of heavy quark limit quenched lattice QCD, we would be reproducing every 
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detail of the results of the Monte Carlo simulations, since that is precisely 
the quantity being measured. 

In our physical picture however, as stated above the transverse gluon 
interactions with the quarks act primarily to renormalize the mass (and in 
the lattice also the Wilson r parameter) in the quark kinetic term and the 
net effect of Coulomb gluon exchange between the light and static quarks 
can be expressed by the potential acting between two infinitely heavy color 
sources. Under these conditions the equation satisfied by qH reduces to: 

<0| (l°E^- t ~ il ■ V + il° E A + m)q H \P*) = (20) 
that together with the Heisenberg equation ^.qu = [H, qu] (in Euclidean 
space) and the relation (0|[i?, quWP*) = —Ebs{^\qh\P*), give rise to the 
eigenvalue equation 

{-id ■ V + mfi + V(r)) X (r) = E BsX (f) (21) 

which is nothing but the Dirac equation for the light quark in the presence of 
the confining external field . This equation corresponds, on the lattice, (with 
the renormalization of the Wilson r parameter also taken into account) to 
an effective lattice Hamiltonian given by the usual Wilson fermion action: 

#3 = ]T{x + (x)(m + 3r)/3x(x) 

X 

i 3 

+ ^E^ + ( X + ^) a k X(x) -X + (x) a k x(x + k)] 

3 

£[X+( X + k) X (x) + X + (x) (3 X (x + k)] 
z k=i 

+ X + (x) F(x) X (x)} (22) 

where x represents a point in the three dimensional lattice of size L , (5 and 
au are just the Dirac matrices, x( x ) is the 4-component wave function, and 
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V(x) is the confining potential determined by Monte Carlo measurements 
of Wilson line correlations of static color sources |IJ . The constituent mass 
m and the Wilson r parameters are free parameters. 

The Hamiltonian #3 defines our new model. From the above discussion 
we realize it represents the closest possible model to QCD consistent with 
the valence quark picture whose validity in the heavy-light meson system 
we want to check. 

As will be shown in Section (3.1) this new model represents a further 
improvement in the prediction of the correct wave functions, that by now 
are, within the errors of the Monte Carlo calculations, essentially fully re- 
produced, indicating the validity of the valence quark model to describe 
heavy-light mesons. Given the necessary assumptions to generate this pic- 
ture from QCD (stated above), the strong coupling nature of the confining 
mechanisms, and the lightness of one of the quarks clearly reflected in the 
necessity of a fully relativistic kinetic energy, the success of the model can 
hardly be expected a priori, and constitutes a strong statement about QCD 
dynamics. 

3.1 Comparison with data 

To find numerically the eigenvectors and eigenvalues of H3 , although we fol- 
lowed in general the same procedure outlined in Section(2.2), some features 
of H% had to be taken into account. For example, due to the non-positivity 
of the spectrum, the inversion of E — H3 was performed with a generalization 
of the conjugate gradient algorithm, the so called minimum residual algo- 
rithm 0|, that takes care of symmetric but non-positive definite matrices 
(one may replace the N x N complex hermitian H3 with a real symmetric 
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2N x 2N version, which is however non-positive-definite). To locate the 
correct region of the spectrum we started in the large mass regime where 
the wavefunctions are well understood and gradually reduced the mass while 
tracking the resulting eigenstates. 

In the following figures we present the results of our new model and 
compare them with the Monte Carlo Data and the predictions of H2. The 
values of the parameters are chosen again to reproduce as well as possible 
the ground state of the system. Following as closely as possible what was 
done in the Monte carlo simulations, (briefly described in section 3), the 
results of H3 presented in the figures, constitute the average of the two 
upper components of the corresponding four-component eigenvectors. 

In FigQ we see that our new model performs as well as H2 for the 
ground state, where there was essentially no room for further improvement. 
We should however note that while the optimum value for r in H3 suffers only 
a small change with respect to the one in H2, the optimum mass becomes 
considerably heavier. 

Once the parameters have been fixed to reproduce as well as possible the 
ground state of the system, we may compare the IP state. Again, as in the 
previous figures for IP wave functions, we divide them by cos and present 
only the radial part. In this case we clearly see the quantitative superiority 
of H3 over the previous models. Near the origin H3 falls much closer to the 
data than Hi- 

We see then that choosing the optimal parameters for the respective 
models, a full Dirac model based on the operator #3 (that, as we have seen 
in section 3 is conceptually as close to lattice QCD as possible within the 
valence quark model), outperforms all the other models and within Monte 
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Figure 7: IS state, L = 12. The Dirac model performs in this case slightly 
better than H2, although there is little room for further improvement in this 
case. 
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Figure 8: IP state, L = 12. The Dirac model performs in this case much 
better than H2. 
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Table 1: 



Model 


E2S — Eis 


Monte Carlo, re = 0.168 


0.31 ± 0.02 


Hi, m = 0.23 


0.381 


H 2 , m = 0.23, r = 0.54 


0.356 


Hz, m = 0.4, r = 0.5 


0.324 



Carlo errors essentially fully reproduces the data. 

We had also available the energies of the IS and 2S states for the Monte 
Carlo data, obtained from the multistate smearing analysis of The 
only meaningful comparison is between energy differences since there is an 
arbitrary choice in deciding the zero energy of the potential V(r). The 
respective energy differences between IS and 2S states are presented in Table 
1. 

Again model H% is in better agreement with the Monte Carlo results 
than the others and, within the errors, reproduces the measured results. 

Model H3 was systematically closer to the data for other values of the 
hopping parameter re. We present in Table 2 the energy splitting for the 
Monte Carlo data corresponding to re = 0.161. This value corresponds to a 
heavier light quark and the optimum values of the parameters correspond- 
ingly change. They are also presented in Table 2. Although the Monte Carlo 
predictions for the various values of the energies change with respect to the 
previous ones, the energy difference essentially remains unchanged. This 
behavior is closely followed by H3 that continues matching the data. Very 
interestingly though, H2 suffers an appreciable modification in the right di- 
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Table 2: 



Model 


-^25 — E±S 


Monte Carlo, re = 0.161 


0.32 ± 0.02 


H 1 ,m = 0.32 


0.385 


H 2 , m = 0.32, r = 0.46 


0.338 


H 3 , m = 0.5, r = 0.45 


0.325 



rection, it's predictions approach the ones of H3 for this heavier case. The 
approach of models 2 and 3 for heavier quarks will be dsicussed in greater 
detail in the following section. 

In the next section we will discuss the nature of the improvement of H% 
with respect to the previous models. 

4 Physical Origin of differences 

In order to fully appreciate the nature of the quantitative improvement given 
by our new model, we will now compare it with the previous ones. 

An obvious difference between the model given by Eq.(^2|) and those 
described by equations (1) and ( |i"3"D is that the former takes into account spin 
effects. The Monte Carlo wavefunctions with which we have tested the model 
were in fact spin-averaged, but H3 contains in principle a full description 
of spin-orbit effects. What follows is a comparison of the models at the 
spin-averaged level. The MonteCarlo wavefunctions obtained in heavy-light 
simulations are typically obtained by averaging the two upper components 
of the light quark propagator on the final time slice. That is why we have 
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performed the same averaging when computing a meson wavefunction from 
the new RQM. 

Expressing the kinetic part of the Hamiltonian H% in momentum space, 
we get, 

1 3 
H 3 kin = 73 E X + (q){M(q) (3 + ]T a k Q k (q)}x(q) (23) 

L q k=l 



with 



3 

M(q) = m + r^(l -cosq k ) (24) 
fe=l 

Qk(q) = sing fc (25) 



Observing Eq.([23|) and Eq.(lS), we realize that a meaningful comparison 
requires expressing the Dirac- Wilson Hamiltonian of Eq.([22]) in a represen- 
tation in which the kinetic energy acquires the form of the kinetic energy 
piece of Eq. (O) . In the continuum this representation exists and is given by 



the well known free Foldy-Wouthuysen (FW) transformation [11]. By this 
we mean a transformation where the Dirac field is rotated by the unitary 
transformation which decouples upper and lower components in the absence 
of interactions. Of course, the full Foldy-Wouthuysen transformation per- 
forms this decoupling including the interaction with the external field order 
by order in the inverse quark mass. However, we wish to avoid a large mass 
expansion for light quarks, and an "all-orders" version of the FW transfor- 
mation is not known explicitly. Nevertheless, the relation between models 
H2 and H3 can still be clarified by a partial FW transformation in which 
upper and lower components are decoupled in the kinetic term only. On the 
lattice the corresponding representation goes along the same lines as in the 
continuum. We then write 



29 



H' 3 = x + e- tb e* b H 3 e- lb e ib x (26) 

where e lS is a unitary (but nonlocal) operator . In momentum space, if we 
choose e lS according to (See [4]) 



(p|e i5 |q) = L 3 <) Piq [cos9 q + ^; ( ^ sin6 c 



(27) 



where Qj(q) is given by Eq.(f25|), 7* are the Dirac gamma matrices, and 



COS 0q EE -L 



sine q EE i= 



\ 



1 + 



1 + 



IQ(q)| 2 



\ 



1 + 



IQ(q)l' 



(28) 
(29) 



A/*(q) 

After this transformation, the kinetic part of H 3 becomes 

^3kiu = ;^£* + (P)/^*(P) 



(30) 



where E p = ^/M 2 (p) + ELl Q*(q), with M (p) and Q*(<0 S iven b Y (0) 
and (|25|), and ^ = e l,s x- 111 Ref , a lattice Foldy-Wouthuysen transfor- 
mation is also being considered. 

So now, both models have the same kinetic part and the difference be- 
tween them becomes completely transparent. Namely, while the model of 
Eq.(li~3|) has (in coordinate space) a potential energy of the form: 



H 2pot = ^+(x)V(x)V(x) (31) 

X 

the potential energy of the model H% becomes after the Foldy-Wouthuysen 



transformation of Eqs (26-p 
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^3 pot = ^E* + (^{EE emz ~ x;eW ~ w ^) 

z x,y p,q 

[cos 6 q + sin 8 q ] [cos 9 P - sin 6 p ]}tf (y) (32) 

\Q\ \f\ 

as can be seen simply by expressing the fields x an d X + m terms of ^ and 
ty + through x = e~ tS ^ and x + = ^ + e lS . 

Comparing Eqs (|3l| ) and (|32] ) and taking into account the definitions of 
cos P and sin p given by Eqs ( |28| ) and ( ]29| ) , we see that p2] ) reduces to 
(31) in the m — > oo limit, in which cos @ p — > 1 and sin p —* and therefore 

^ pot m -^° ^E^IeE^ - ^ 



= E^ + ( f ) y ( f )^( f ) ( 33 ) 

It is worth looking at the above limit in more detail. Expanding the 

product between brackets in Eq(|32|) and remembering that = — ia % \ 
we obtain 



^ pot = ^E^ + (^{EE e *" xle ^ v ^) 

z %,y p,q 

[Fi (p, q) + F 2 (p, q) + F 3 (p, q)}}^(v) (34) 



where 



Fi ip, 0) = cos 9 q cos 9 P + |Q^||p|pj| sin ®q sin ®p ( 35 ) 
= ial3 | Q'/q) |^(p)| sine q sine P ( 36 ) 

(p, q) = i % |q^| sin ®q cos p ~ |p^p)| sin 0p cos 0q ^ 
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To interpret these terms it is convenient to consider their continuum 
limit. In this limit, Eq(|2^) and ( p9| ) become: 

l|q| 2 



COs6 q EE -L 



sin q EE — 



\ 



1 + 



1 



m— >large 



1 



m- 



V2\ 



1 



1 



♦large 1 |q| 



'l i |q| 2 2 m 

and the interactions corresponding to the three terms above become: 



(38) 
(39) 



ix,y Jp,q 
|2 1 |„|2 



[i-^^V- S 1! T + 7 1 fW ( 4 °) 



1 !s!_. _ 1 IpI , iq-Pn 

m 2 8 m 2 4 m 2 
^ + (x)F(x)*(x) (41) 

+ -^3 J^ + (x)V 2 V(x)*(x) (42) 

Term (f4l|) represents the electrostatic energy of a point-like particle and is 
the one present in models Hi and H2. More interestingly, term (42) corre- 
sponds to exactly the Darwin term. It arises because of Zitterbewegung, as 
can be seen from the smearing of the potential (see Eqs(|34|) and ( |35| ) ). 



H 3 f ot = [$> + (z){[ [ e^-VjW-ftVix) 

J z J x,y Jp,q 

[-^7^%W (43) 
L |q||p|4 m 2 U yy ' K ' 

= [y + (z){[ [ e^-^e^-^Vix) 

J z J x,y Jp,q 

[«««»( o < 44 > 

= -i//^^^ (45) 
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Where we have used the identity a lJ = e^ 



a 







a k 



, where the a k are 



the Pauli matrices. Clearly Eq(45) represents the spin-orbit interaction. 
Finally we have: 



H 3 % = l^ + {z){! [ e^-^e^-^Vix 



x,y Jp,q 



[7^-7^]}*^ (46) 
2m 2m 

= 2^X* +(f){ ^ F(f)7 ' } * ( ^ (4?) 

representing interactions between upper and lower components of the Dirac 
spinor. In a large mass expansion (the usual FW transformation) this term 
is removed by a unitary rotation at order 1/m 2 . 

It is perhaps worth mentioning that ignoring this term completely (clearly 
valid for large masses only!) but not making the large mass expansion in ( fi"0| , 



43) and spin-averaging the resulting Hamiltonian, yields a modified poten- 
tial model which we have studied and which yields wavefunctions very close 
to the full Dirac formalism. This approximation is not very well motivated 
however, as it seems to involve a rather inconsistent treatment in terms of 
a 1/m development. 

In retrospect we realize that potential model descriptions based on H\ 
and H2 are somewhat inconsistent since, as we have just seen, they effectively 
take the limit m — > 00 in the potential part while keeping a finite mass in 
the kinetic part (as first shown in Ref Q , the full relativistic kinetic energy is 
essential in reproducing the data). The reasonably good agreement between 
Hi, H2 and the Monte Carlo data, together with the inconsistency pointed 
out above, deserves some comments. The validity of H3, as clearly stated 
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in section 3, is based on the assumption that, even though the light quark 
moves fast enough for relativistic effects to be important, the time scales 
over which, the string connecting the quark to the static source responds to 
changes in the light quark position, are small compared with the time scales 
relevant for the light quark motion. This makes the interaction between the 
light and the static quark well described by the energy which would obtain 
if the light quark were held fixed. Models Hi and H2, effectively taking the 
limit m — > 00 in the potential part and keeping a finite mass in the kinetic 
part, are simply making the further assumption that the confining potential 
is essentially constant over regions of size of the order of the light quark 
Compton wave-length. To see this implication we just have to remember that 
a Dirac particle does not move along a straight line with constant velocity 
but instead carries out an oscillatory motion (Zitterbewegung) with the 
speed of light (see [11, [l3|]) centered on a point which does move uniformly. 
This oscillatory motion is of the order of the Compton wavelength of the 
particle. As our light quark moves though the confining potential, its color 
charge explores then the field over a region of the order of its Compton 
wavelength and this explains the appearance of the Darwin term and all 
higher order terms familiar from the F-W transformation. However if over 
regions of the order of the Compton wavelength the field is slowly varying, 
it may be reasonable to ignore the smearing effects (formally higher order 
in 1/m) while maintaining the relativistic kinematics in the kinetic term. 
This seems to be the case in our situation in which, as can be seen from 
the reasonable success of models H% and H2, taking the m —* 00 limit in 
the potential part seems not to be a very bad thing to do (for example, at 
larger light quark masses than those studied here, the agreement between 
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the wavefunctions generated from Hi and #3 is closer). However were we 
going to do the same in the kinetic part, we would get a non-relativistic 
model that does a very bad job at reproducing the wave functions Q. 

In any case, the most important lesson that we learn from the discussion 
above is that the differences that we saw in the previous section between the 
wave functions of model 2 and those of model 3 are, as we have just seen, the 
result of well known effects that arise when one combines quantum mechanics 
and relativity, which model 3 captures (to the extent that the Dirac equation 
captures them), but are ignored in models 1 and 2. These effects are to our 
knowledge visible for the first time in the context of quantitatively measured 
(in quenched lattice QCD) strong interaction wavefunctions. 

5 Conclusions 

We have presented the results for the IS and IP wave functions and energy 
differences between the IS and 2S states of a fully relativistic lattice model 
of heavy-light mesons. These results were compared with Monte Carlo mea- 
surements of the corresponding quantities and with previous models. The 
results of the comparison validated the valence quark model as a good rep- 
resentation of heavy- light mesons, at least for the lattice sizes tested so far. 
In particular our fully relativistic model proved quantitatively as well as 
qualitatively superior to previous models. The quantitative improvement 
represented by our model arose simply by comparison with the data. The 
qualitative one came not only from the relative transparency of the ap- 
proximations being done, clearly stated in the derivation of the model; a 
comparison of the physical content of the different models revealed that the 
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previous ones were somewhat inconsistent in their relative treatment of the 
potential and kinetic terms. It is precisely this comparison that allows a 
physical interpretation of the quantitative improvements of the fully rela- 
tivistic models. As it turned out they can be thought of as due to Darwin 
and higher order effects (in the language of a Foldy-Wouthuysen treatment) 
arising from the quantum-relativistic derealization of the light quark due 
to Zitterbewegung. It is remarkable that the Monte Carlo simulations of 
Ref [pj are now accurate enough to capture this phenomenon. 

We expect to be able to extend the above results to much larger lattices. 
We are currently generalizing this work to treat mesons with two finite 
mass quarks. If the fully relativistic model continues to be as quantitatively 
accurate as the results obtained here suggest it may turn out to be a very 
useful tool in the study of the spectrum and static properties of charmonium 
and charmed and B-mesons. 
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